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Abstract. We study the problem of finding an optimal radiotherapy treat- 
ment plan. A time-dependent Boltzmann particle transport model is used to 
model the interaction between radiative particles with tissue. This model al- 
lows for the modeling of inhomogeneities in the body and allows for anisotropic 
sources modeling distributed radiation — as in brachytherapy — and external 
beam sources — as in teletherapy. We study two optimization problems: mini- 
mizing the deviation from a spatially-dependent prescribed dose through a qua- 
dratic tracking functional; and minimizing the survival of tumor cells through 
the use of the linear-quadratic model of radiobiological cell response. For each 
problem, we derive the optimality systems. In order to solve the state and ad- 
joint equations, we use the minimum entropy approximation; the advantages 
of this method are discussed. Numerical results are then presented. 



1. Introduction 

Radiotherapy is one of the main tools currently in use for the treatment of can- 
cer. Radiation is deposited in the tissue with the aim of damaging tumor cells and 
disrupting their ability to reproduce. In this paper, we are interested in the treat- 
ment planning problem. We wish to determine a method of delivering a sufficient 
level of radiative energy to the tumor cells to ensure cell death. This is balanced 
by our our desire to minimize damage to healthy tissue, especially specific critit- 
cal structures — which we call regions at risk — that should be damaged as little as 
possible. This problem can be divided into two phases. First, one must deter- 
mine an appropriate external source of radiation that optimally balances these two 
competing goals. Secondly, one must devise a means to administer this external 
source through (in the case of teletherapy) the selection of beam angles and sim- 
ilar parameters. We focus on the first of these steps. Once the optimal source is 
determined, the means of creating such a source would still need to be obtained. 
Often, the optimal dose will be determined manually and iteratively by an expe- 
rienced dosimetrist; this dose consists of selecting several fixed beam angles and 
strengths. However, newer methods such as Intensity Modulated Radiation Ther- 
apy allow for dynamic sources; such a method requires more automated algorithms 
for determining the treatment scheme |13j . 

Quantifying what the "best" treatment plan is difficult and still unsettled (see 
[13], [H] for instance). Often a heuristic approach is used by technicians. Here, 
we will study two different means of measuring the effectiveness of a given source. 
First, we use a general measure: we look to minimize the deviation from a desired 
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distribution of radiation, similar to the objective functional used in [13j . [2].|10j. 
and [H]- The other method we focus on involves using radiobiological parameters 
describing cells' susceptibility to damage from radiation. The cells' response to 
the delivered dose is modeled by the frequently-used linear-quadratic model for cell 
death. A significant difficulty in using this model is that the parameters for cell 
damage are dependent on several factors leading to significant levels of interpatient 
heterogeneity [H] , [IH] , [H] ) [S] ■ Furthermore, the parameters can change over the 
course of treatment plans as patients exhibit changing levels of radioscnsitivity. In 
particular, [5^ studied the optimal treatment planning problem and the dependence 
of the optimal dose on different parameters, especially in the face of changing 
radioscnsitivity. In [,14 , a method for treatment planning that accounts for the 
uncertainty in an individual's specific linear-quadratic parameters was developed. 

As described in [10) , the Fermi-Eyges theory of radiation is used for most dose 
calculation algorithms; this theory fails to incorporate properly physical interactions 
in areas where inhomogeneities such as void-like regions occur. Here, we will make 
use of a Boltzmann transport model which describes the scattering and absorb- 
ing effects of the tissue while also allowing for inhomogeneities. Our optimization 
method here differs in two significant ways from that in |TU1. First, our model is 
time-dependent; we will allow for time-dependent, anisotropic sources which are 
external (teletherapy) as well as sources distributed in the tissue (brachytherapy). 
Second, as opposed to the spherical harmonics approximation used there, we will 
make use of a minimum entropy approximation to the Boltzmann transport model 
which has the benefit of maintaining several desirable physical properties from the 
original model. 

In Section [2j we will introduce the Boltzmann transport model describing the 
interactions of the radiation with the tissue and the two objective functionals which 
we use to measure the effectiveness of a given dose. In Section [3] we establish 
optimality conditions for both objective functionals through the adjoint equations 
of the Boltzmann equation. From there, we introduce the minimum entropy (Mi) 
approximation and briefly discuss the advantages of this approximation. Numerical 
results in two dimensions are presented for several test cases in Section [sj 



2.1. Transport Equation. Let Z C M'^ be an open, convex, bounded domain with 
smooth boundary and outward normal n{x), and T > be given. The direction of a 
particle's movement is f2 € S'^ where S"^ is the 3-dimensional unit sphere. Then we 
model particle transport with sources : [0,r] X X 5*2 ^. M by the following 
equation 



V'(t,a;,Sl) ^ qi{t,x,n), \f{t,x,n) e F" := {{t,x,n) : t G [0,T],x e dZ,n{x)-n < 0} 



2. MATHEMATICAL MODEL 



(1) 

tpt{t, x,n) + n- \7^^p{t, X, n) + at{x)ij{t, x, n) 



+ q{t,x,n) 



= as{x) [ s{x,n-n'yj{t,x,n')dn' 




and initial condition 



(3) 



ip{o,x,n) = 0, v(x,f7) ezx s^. 
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It should be noted that, ([T]) does not explicitly include the energy dependence of 
the transport of particles. However, as detailed in [9], the Boltzmann Continuous 
Slowing Down model, which does include energy dependence, can be transformed 
into a problem such as ([T]) where the initial- value condition intuitively describes the 
physically reasonable condition that no particles are of "infinite" energy. In general, 
the term q denotes a distributed radiative source which would be appropriate for 
applications in brachytherapy; meanwhile, qi seems a natural means of modeling 
the effects of external beams in teletherapy applications. However, as described in 
Section [4j the role of the control on the incoming boundary is difficult to model in 
the context of the Mi approximation (along with similar moment approximations) 
to ([T]). With this in mind, throughout this paper we will assume qi = 0; that is, we 
will replace ^ with the incoming boundary condition 

(4) ij{t,x,n) = o, y{t,x,n) e r-. 

We will address in Sectiorj4] a means of modeling the case of teletherapy using 
the distributed control. The quantity tlj{t,x,^) cos 6 dAdil is interpreted as the 
number of particles at time t passing through dA at x in the direction dVl near f2 
where 9 is the angle between and dA. The total cross-section at{x) is the sum 
of the absorption cross-section aa{x) and scattering cross-section as{x). Finally, 
the scattering kernel — intuitively describing the probability of a particle changing 
direction from Q' to fl at x € Z after a scattering event — is normalized; that is, it 
satisfies for all x, 

In this paper, we will usethe simplified Henyey-Greenstein kernel (see [J for more) 

( \ = 1 - 9{^? 

^hgKx, t]) . ^^^^ ^ ^^^^^ _ 25(^)77)3/2 

where g is the average cosine of the scattering angle. Following [TO], we define the 
operators 

iA^){t, X, n) = -Vt(i, x,n)~n- v^^{t, x, n) 
{Ki;){t,x,n) ^ as{x) [ s{x,n' ■ ^})^p{t.x,n')dn' 
{i:^P)it,x,n) = at{x)ij{t,x,n) 



and make the following general assumptions: 

(H-1) (Tt,CT„S>0, 

(H-2) cTt,as€ L°°{Z), 

(H-3) CTtix) -'cr,{x)>a> 0, Vx G Z. 



We note that (H-I) and (H-2) are reasonable physically and that as long as 

is 



the medium under consideration is absorbing, the coercivity hypothesis (H-3) 
satisfied. These assumptions imply that the operators 

K,S : L2([0,T] X Z X 5*2) ^ i^([0,T] x Z x S^) 

are bounded and linear operators. Finally, we note that if we define 

D{A) = {t/- e L^{[0,T] X Z X S"^) : il)t,n ■ e L^i[0,T] x Z x S"^)}, 
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then 

A : D{A) L2([0,r] x Z x S^) 

is a well-defined linear operator. We also let D{A) denote the set of functions of 
D{A) which also satisfy Q and (|3|, for qi = 0. Then we can conclude — for instance, 
by Ch. XXI, 2, Theorem 3 of [5 the following holds. 

Theorem 2.1. Suppose that \(H-l)f(H-3)\ hold and that qi e L'^i[0,T] x Z x S'^) is 
nonnegative. Then ^ has a unique solution ip C D{A). 

We define the control-to-state mapping 

8 : L2([0,T] xZxS^)-f D{A) 

as the map which maps a source q to the corresponding solution ^ to 

-All; + J:tp = Ki/j + q 

subject to Q and ([s]). By the above, this is a well-defined, linear operator. As in 
|10) . it is also bounded. 

Proposition 2.1. Let (•,-)2 denote the standard inner product on L^([0,r] x Z x 
5^). Then it holds that, for ^ satisfying Q) and 



12 

and thus £ is a bounded operator. 

Proof. We note that, by an argument analogous to the proof of Lemma 3.1 of [10], 

(S^-K^,V')2>a||^||i 
Additionally, by integration by parts and the initial and boundary conditions, 



-AV', -0)2= / {'ijjtip)dtdxdn + / {n ■ V xi^)ipdtdxdn 

= / l-'tij^{T,x,n)dxdn+ / / div( / lni:^dn)dxdt 

JzxS^"^ J[0,T]Jz ^Js^'^ ' 

^'iP^{T,x,n)dxdn+ [ [ -{n-n{x))+^/j^dxdndt>o. 



ZxS^ 2 J[0,T] JdZxS^ 2 

Combining the two inequalities gives the desired result. □ 

2.2. Objective Functionals. We consider two objective functionals. In each case, 
we neglect the incoming boundary source term by fixing qi = for reasons which 
are discussed in|4] Also, we are primarily interested only in the total dose deposited 
as a function of x, which we denote by the operator D : L'^{[0,T]x Z x S^) — > L'^{Z) 
which is defined as 



:= / i}}{t,x,n)dtd9.. 

J[0,T\xS^ 

The first functional is a quadratic tracking which measures the deviation from a 
prescribed dose D{x) 

(5) Jr(V','7):= / ci{Di,-Dfdx+ f C2D{qfdx. 
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Here ci,C2 > are spatially dependent weighting functions. One such weighting 
function, used in |10j . for ci is 

Cl = CtXZt + CRXZa + CNXZm : 

where ct, cr, cn are constants and Zt, ^_r, Zn are, respectively, the portions of Z 
corresponding to the tumor, a region at risk, and normal tissue with Z = Zt U 
Zn U Zn. This functional differs slightly from that studied in jlO] in that we are 
concerned with the dose deposited but the source is allowed to be angle-dependent. 
In that work, the objective functional studied for anisotropic controls also assumed 
a desired anisotropic distribution was to be tracked, as opposed to simply a total 
dose. Though this objective functional does not attempt to model cell response 
to the dose, it does have the benefits of not requiring as inputs the radiobiological 
parameters needed for such a model. 

The second objective functional involves modeling cell death as a function of 
D{ip). We will use the linear-quadratic model neglecting cell growth and 

repair rates. In a reference volume for a given dose D, the fraction of surviving 
cells of a single type is given by 

(6) SF = exp(-ai:i - I3D^) 

where a,/3 are parameters dependent the type of cell being irradiated. This model 
has been used extensively for determining dose strategies (for instance, [il lTClfT^ITSl 
lllj and, in a modified form, [S]). The surviving fraction is often used to determine 
the tumor control probability for a given region: 



TCP = expj ( - pexp{-a,D, - P.Df^ 



where p is the density of tumor cells as a function of space. We consider a continuous 
measure of the number of surviving cells of the various cell types in Z. We consider 
a single type of tumor cell and assign this cell type the index 0. The remaining 
cell types are given index i = 1,...,N. Then, for cell densities pi{x) > with 

J2i=oPii^) = 1' define 

(7) Jsf(^,9) flo / poexp[- aoD^: - i3o{D^/j 



Z 

N „ 

J2a, y^ft(l-exp[-a,DV'-APV')']) 



i=l 
C2 

2 



[Dqf 



z 



where we assign constant weights for each cell type tti > and the control C2 > 0. 

For either cost functional, we set U{x) > to be the maximum allowed source 
at X e Z. Let L'^^([0,T] x Z x S"^) denote the set of admissible controls-that is, 
those controls satisfying q > and 

q{t,x,n)dn < U{x) 

almost everywhere. Note that this subset of i^([0, T]x Z x S'^) is nonempty, closed, 
and convex. 
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3. OpTIMALITY CONDITIONS 

We are, then, interested in two separate optimal control problems: the mini- 
mization of tracking error 

{Vt) min Jrii'.q), s.t. £{q) = ip. 

ipeD{A),qeLlJO,T]xZxS'^) 

and the minimization of tumor cell surviving fraction balanced against normal cell 
death 

(Vse) min JsF{tp,q), s.t. £{q) ^ tp. 

jpeD(A),qeLl^{[0,T]xZxS^) 

In order to obtain optimality conditions for Vt and Vs f , we follow the example of 
|17) . One of the advantages of using Jt as the objective functional is the following: 



Theorem 3.1. Given \(H-l)'l{(H-3)\ and V' G D{A) with D = D{ip) and ci,C2 > 0, 
then Vt has a unique solution q* G L'^^{[0,T] x Z x S"^). 

Proof. We note that the reduced objective functional jT{q) — Jt{£{(1), q) is a con- 
vex functional over the closed, convex set L^^([0,r] x Z x S'^) of the Hilbert space 
L'^{[0,T] X Z X S"^) and that £ is linear and bounded. By Theorem 2.16 of [II], 
there exists a unique minimizer for problem Vt- D 

We now turn to the existence of an optimal control for Vsf- Because the asso- 
ciated reduced objective functional is nonconvex, we can not assume existence of a 
unique optimal control. 



Theorem 3.2. Given (H-1 )^H-3) and Ci,a > there exists at least one optimal 
^^^i^^l ^* ^ r2 /rn T^l V 7 V Q'^\ f^^ 'D„„ 



control q* e Lla{[0,T] X Z X S^Yfor Vsf 

Proof. We note that Jsf is convex with respect to q and that the set of admissible 
controls is closed, convex, and bounded. This, coupled with the boundedness of the 
operator £ allows us to conclude that the infimum of Jsf over ^^^([0, T]x Z x S"^) 
is finite and that any weakly convergent minimizing sequence of controls (?„ will 
converge to an admissible control q* . By the linearity of A, E, K, we conclude that 
the corresponding solutions converge strongly to a solution corresponding to the 
q* . (see Theorem 5.7 of [H]). The convexity of Jsf as a functional of q allows us 
to conclude this control is optimal, as in Theorem 4.15 of [T7]. □ 

We next note the following from [S]: 
Proposition 3.1. Under \(H-l)}(H-3)\ for any r £ L'^{[0,T] x Z x S'^), 



AA + EA - KA = r 

A = on r+ 

A(r, X, n) = 0, v(x, n) e z xs'^. 

has a unique solution A G L'^{[0,T] xZx S"^), and the solution operator £* is linear, 
bounded, and the adjoint to £. 

Now, we are ready to turn to the optimality conditions for Vt, which is similar 
to that found in [TU] . 
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Theorem 3.3. For a desired distribution D G L^iZ) and under the assumptions 
(H-1 %(H-3) q* £ ^add^'-^] Z X S'^) is the unique optimal solution to Vt if and 



only if 
(8) 



q*(t,x,n)= proj ( ^--\*{t,x,n)) 

[0,(7(2;)] ^ C-)(X) J 



Clix) 



where A* = E* {c\(D{^*) — D)), -j/)* = £{q*), and proj5(-) denotes projection onto 
the closed set S. 

Proof. The gradient of jr is ci£* {D£{q*) — D)) + C2Dq and so we conclude, using 
the same arguments of Theorem 2.25, Lemma 2.26, and Theorem 2.27 of [TJ- Q 

This means that the optimahty system for Vt is 

'-AV'*+I]V* =KV'*+g*, 
(9) iAA*+EA* ='K\* +ciD{£{q*)~D), 

X =Proj[o,t/(x)](-i-^*) 

with V* satisfying Q and ([s]). A* = on r+, and A(r, x, 17) = on Z x 5*2. 

The necessary condition for the optimahty of a source for Vsf is of a similar 
form. 

Theorem 3.4. For a given distribution of cells pi, i — 0,1, ... N and the associated 
linear- quadratic cell death parameters Ui, /3i > 0, a locally optimal q* G L^dii^' ^] ^ 
Z X S'^) for JsF must be such that 



q*{t,x,^)^ proj ( ]—-X*{t,x,n)) 



(10) , 

[0,(7(2:)] ^ C2(a; 

where A* — £*{r*) for 

r* = [( - ao - 2PoDr) (coPo exp[-aoi?i^* - ^0^)^]) 

- E [( - ~ 2/3.^^-*) {c^P^ expl^a^Or - l3,{Drf]) 
1=1 

and ip* = i^^(9*)- 

Proof. The reduced cost functional jsriq) — JsF{£{q),q) is Frechet differentiable 
due to the differentiability of the integrands in Jsf and the linearity and continuity 
of £. For simplicity, we assume po = 1 (and thus pi = for i > 1) first. We then 
get that for a locally optimal q* and any other admissible control q, 

jsF{q*){q-q*)>o 

where, by the linearity of the adjoint solution operator £* and the chain rule. 



fsFillil -q*) = co {-ao~ 2poD£{q*)) (cqPo exphao^^^ (<Z*) - MD£{q*)f])£{q 
J z 

+ j q*iq-q*) 

£*ii^ao - 2l3oD£{q*)){coPocxp[-aoD£iq*) ^ l3o{D£{q*))^]){q - 

+ C2 I q*{q-q*)- 



RICHARD BARNARD* , MARTIN FRANK, MICHAEL HERTY R.WTH AACHEN UNIVERSITY, AACHEN GERMANY 



This leads to the minimum principle (by an argument analogous to that found on 
pages 68-70 of JJTj), which states that 

min [(A* + C2(7*)(7] 

qeLl^([a,T]xZxs^) 

is attained almost everywhere at q* which leads us to conclude that q* satisfies 
(10 1. Now if we drop the assumption that = 1, and allow for the presence of 
other species, the argument follows through in a nearly identical fashion. □ 

This gives the necessary optimality system 

(11) <AA*+SA* =KA*+r*, 

U* = Proj[o,c/(x)](-s^A*) 

where r* is as in Theorem [331 



4. Minimum entropy closure 

For our purposes, solving ([I]) in full is unnecessary. We note that both Jt and 
JsF do not need the exact photon densities as functions of the angular variable; it 
would be enough to determine the zeroth moment, which is the total energy at a 
given {t,x): 

i/,(o)(t,x) = / ii;{t,x,n)dn. 



In light of this, we average ([!]) over the angular variable and obtain 

(12) [t. x) + V, • ^(1) (i, x) + a* (x)^(°) (t, x) = a,(a;) (i, x) + (t, x) 

where 

i!^^\t,x) = I nij{t,x,n)dn 



is the first moment, or the flux vector and g^"-* is the average of q over all directions. 
By multiplying ([I]) by O and taking the average over all directions, we get 

(13) dtij^^'>{t,x)+W, ■ i^^^\t,x) + at{x)^^'-Ht,x) = a,(x)#(i)(t,x) +g(i)(t,x) 

where 

ip(^){t,x)^ {n(g}n)ij(t,x,n)dn 



is the second moment, or the pressure tensor of the radiation field, and g*^^' is the 
first moment of the control. We note that if we continued this process of multiplying 
([T]) by monomials of fl and integrating over S'^, we would always have an equation 
relating the n*'* moment to the (n + l)*^ moment. If we want a closed system, then, 
we must select an approximation to the {n + lY^ moment 

Here, we consider the minimum entropy closure (Mi) to approximate ■0^^). This 
involves searching for / which minimizes 



where 



h*„{I) = 2fcz^2(nlogn- (n + 1) log(n + 1)), with n = — 
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and k is the Boltzmann constant, h is the Planck constant, and v is the frequency 
of the radiation. We constrain the entropy minimization by requiring, naturally, 
that 

/ = and / nidQ ^ ^//^\ 

in [5] , the entropy minimizer was obtained explicitly by setting 

V-^^) - D{f)E 

where the relative flux is / = the Eddington tensor D, is 

^W) 2 ^ 2 |/|2 

and the Eddington factor x is 

x(/) = ^(5-2v/43W)- 

This choice of the Eddington factor preserves several reasonable physical properties. 
Importantly, the flux is hmited; that is, |/| < 1. This preserves the finite bound 
the speed that information travels through the system. Also, the distribution is 
nonnegative, which we would like to ensure. We note that neither property is sat- 
isfied for several common closure choices which use a diffusion approximation, such 
as the spherical harmonics method(also known as Pn methods) [3]. In [7], the Mi 
approximation was seen via numerical tests to provide a high quality approximation 
when compared with Monte Carlo-based methods. 

Despite these advantages, the AIi does have several drawbacks. Due to the non- 
linearity of the Ml approximation (as opposed to the Pn approximation (10)). we 
can not expect that the systems obtained by applying the Mi approximation to 
the optimality systems will match the optimality systems obtained by optimizing 
the Ml approximation to ([!]) . Also, as mentioned briefly above, the use of moment 
models complicates the incorporation of boundary conditions. The boundary con- 
dition Q describes only the incoming particles, whereas boundary conditions for 



(12 1 and Jl3^ require us to prescribe values for the full moments. In light of this. 



we use U{x) to approximate boundary control by setting it as 

if doz{x) < e 
I S otherwise, 

where Qmax > is a total maximum source level, e,d > are sufficiently small 
parameters, and dgz denotes the usual distance from the boundary of Z. We note 
that the first two equations in both ^ and ( 11 ) are very similar and we thus use the 



Ml closure to solve the state and adjoint equations. That is, we take an optimize- 
then-discretize approach to Vt and Vsf as opposed to a discretize-then-optimize 
approach. 



5. 2-D Numerical Tests 

In this section, we study the 2-dimensional region Z := [—1,1] x [—1,1]. We 
define a void-like region Zy ■= {0.8 < |a;| < 0.9} On Zy, the scattering and 
absorbtion cross-sections as, da have values that differ from those elsewhere; these 
are summarized in Table [T] for the specific values. 
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Zv 


Z\Zv 




0.001 


0.05 


OS 


0.01 


0.5 



Table 1. Material parameters for Vsf 







Zr 


Zn 


a 


0.52 


0.170 


0.170 


fi 


0.171 


0.0078 


0.0078 



Table 2. CcU-rcsponse parameters for Vsf 







Zr 


Basic Target 


[-0.25,0.25] X [-0.25,0.25] 


[.254,0.379] X [-0.125,0.125] 


Intermediate Target 


([-0.25,0.25] X [-0.25,0]) 
U([-0.25,0] X [0,0.25]) 


[0.04,0.25] X [0.04,0.25] 


Complex Target 


([-0.25,0.25] X [-0.25,-0.125]) 
U([-0.25,0.25] X [0.125,0.25]) 
U([0.04,0.25] X [-0.125,0.125]) 


[-0.25,-0.04] X [-0.121,0.121] 



Table 3. Locations of Tumor and Risk Region 



The final parameter of the physical medium to be defined is the scattering phase 
function. For simplicity in our test problem, we define g = 0.85. Other, more 
involved, kernels might be used in this range for g to agree with Mie scattering 
theory as in |T]. 

For testing, we use three tumor/risk region regions similar to those in . Specif- 
ically, we define the regions in Table [3] and are shown in Figure [ij the void region 
is shown in black and the tumor and risk regions are traced in white. In the basic 
target case, seen in Figure [l (a) the tumor region is a box, as is the risk region. The 



second, intermediate target case, seen in Figure 1(b) involves an L-shaped tumor 
around a box-shaped risk region. Finally, the complex target case in Figure |l(c) 



involves a C-shaped tumor around a risk region. We will solve both Vt and Vsf 
for each geometry seen in Figure [l] For Vt for each example, we set D = T on Zt 
and elsewhere, corresponding to an average (over time) dose of 1, and 

Cl = CtXZt + C-RXZn + CnXZm 

with CT = 25, cr = 150, cn = 1. For Vsf, we set po = xzt, Pi = XZr, P2 = 
XZk, with weights oq = 500, ai = 2000, 02 = 1. The LQ parameters are shown 
below; the tumor region uses parameters for Lewis Lung tumor cells and the risk 
region and normal region use parameters for V79 (normal Chinese hamster lung 
tissue); both sets of parameters are taken from [T^. We use a simple gradient 
descent algorithm to determine an optimal control for each of the three targets 
and for each objective functional. The Mi approximations to both the state and 
adjoint equations are solved with a finite volume solver using approximately 10,000 
cells with final time T = 5. The parameters S and e in the definition of U (x) 
are chosen as min{Ax, A?/}(10^'*)i and min{Aa;, Ay}, respectively, where Ax, Ay 
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(c) Complex Target 

Figure 1. Example cases 



are the mesh sizes. We consider the algorithm as converged if both the difference 
between consecutive iterates and ||j'(g)||oo are smaller than fO~^. 

Figure [2] shows the optimal boundary source term for both Vt and Vs f ■ The 
vectors shown on the boundary ar e the time -inte grate d values of g^^-* normalized 



and then scaled by ' . In Figures 2(a) 2(c) and 2(e) (corresponding to Vt), the 



isolines are spaced at 5% intervals of the maximum of the desired dose (here, 5). 
In the intermediate and tracking cases, we see that relatively low dose levels are 
attained, primarily due to the high penalty to any dose deposited in the risk region. 
In Figures [2(b)[ |2(d)[ and 2(f)[ corresponding to Vsf), the isolines are spaced at 
intervals of 10% of cells killed. Here a high proportion of the tumor cells are killed 
(in each case 80% — 90%) while in the Intermediate and Basic cases, the tumor has 
at least 60% survival; in the Complex case, the risk region has 50% — 60% survival. 

The dose deposited in Vt changes significantly when we alter the relative weights 
ct and c_r. In Figure ??, we see the results for solving Vt, and Vsf where we set 
cu = 50 and oi = 1000 (all other parameters are unchanged). In both the basic 
and intermediate cases, the dose delivered to the tumor is significantly higher while 
also remaining largely concentrated on the tumor. This is slightly less true in the 
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(a) Tracking, Basic Target (b) Surviving Fraction, Basic Target 




(e) Tracking, Complex Target (f) Surviving Fraction, Complex Target 



Figure 2. Optimal boundary source for Vt, in (a)|(c)|(e) and 
VsF in |(b)|(d)|(f)| 
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complex case. In the figures for VsFi there is no general change in the pattern 
of cell death, however the risk region in the first two cases has 50% cell survival 
whereas the complex case has approximately 40% — 50% cell survival. However, 
this lack of change in pattern can partially be attributed to the tumor cells being 
more susceptible to the radiation dose. 

We conclude with a final set of numerical examples which restrict the location 
of the source by altering the definition of U{x). Here we require that q < 5 on 
one side of the boundary. For the basic and intermediate case, we require that the 
external source not come from the left side of Z. For the complex case, we disallow 
sources on the right side (as the optimal source is nearly zero on the right side in 
the complex case for Vt)- Figure |4] shows the optimal solution for both problems, 
using the same penalization parameters used in Figure [2j The optimal dose for Vt 
is significantly worse, with the tumor in the intermediate and complex cases getting 
a dose below 3. However, the tumor cells have a survival of 10% or less for each 
case and the risk region has a survival rate of 60% or higher in each case. 

Clearly, the choice of relative weights ct, cr and affects the optimal solution. 
However, we have currently no systematic a priori method of selecting the relative 
weights in either problem. Such a method should take into account the relative 
importance of the different tissues and the susceptibilities of the various cell types, 
leading to a possibly nonintuitive scheme. Without such a method, several opti- 
mization runs may be required to find an appropriate treatment strategy. We also 
note that the optimization does not consider the feasibility of the optimal source. 
In fact, the optimal boundary control in all of the numerical results are spread out 
along the boundary, as opposed to being focused in a beam configuration. Further 
modeling would be required in order for the cost function to penalize non-beam like 
configurations. 
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Figure 3. Optimal results with lower penalty to dose in Zji 
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Figure 4. Optimal results with one edge blocked 
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